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I.  INTRODUCTION 


A  detailed  understanding  of  the  phenomena  occurring  within  a  gun 
tube  during  firing  would  aid  in  the  development  of  modern  high  muzzle 
velocity  gun  systems.  In  particular,  the  reliable  prediction  of  flow 
properties  in  the  gun  tube  would  aid  in  the  understanding  and  prevention 
of  both  gun  barrel  erosion  and  catastrophic  gun  failure.  The  ballistic 
cycle  which  we  wish  to  simulate  begins  with  the  ignition  of  the  solid 
propellant  particles  and  ends  with  the  projectile  exiting  the  gun  tube 
and  the  discharge  of  hot  gases  and  any  unburned  propellant.  The  igni¬ 
tion  of  the  propellant  particles  is  normally  initiated  by  a  hot  gas- 
particle  flow  from  a  primer  consisting  of  an  igniter  and  a  black  powder 
charge.  The  hot  primer  gas  flows  into  the  propellant  bed  and  heats  the 
propellant  particles  via  conduction,  radiation  and  convection.  As  the 
propellant  burns  gas  pressure  and  temperature  rise,  thereby  promoting 
the  flame  spread  through  the  propellant  bed.  Pressure  waves  propagate 
through  the  chamber,  and  after  the  pressure  becomes  large  enough  to 
overcome  frictional  resistance  to  motion  the  projectile  will  accelerate 
down  the  tube.  At  some  point  in  the  cycle  the  propellant  bed  will  flu¬ 
idize,  and  interaction  among  the  particles  will  diminish.  As  the 
projectile  moves  down  the  tube  several  processes  may  occur  which  are  not 
considered  in  the  present  analysis;  these  are  projectile  instability  and 
tube  vibration  and  leakage  of  propellant  gases  past  the  projectile. 

The  former  two  phenomena  would  not  normally  have  a  significant  influence 
on  the  two-phase  flow  in  the  tube.  While  leakage  of  propellant  gases 
may  influence  the  flow  especially  near  the  projectile-tube  wall  junc¬ 
ture,  its  consideration  has  been  deferred  until  a  later  stage  of  the 
program  development. 

Until  recently,  most  other  efforts  in  modeling  the  flow  phenomena 
in  guns  have  been  limited  to  quasi-one-dimensional ,  inviscid  two-phase 
flow  analyses.  *  ,  *  Portions  of  the  ballistic  cycle  can 
now  be  modeled  with  two-dimensional  analyses  by  assuming  cyl indrically 
symmetric  flow.  Gough5  has  developed  an  inviscid,  two-dimensional 

ZT  P.S.  Gough ,  F.J.  Z warts,  " Modeling  Heterogeneous  Two-Phase  Reacting 
Flow,"  AIAA  Journal ,  Vol.  17,  Ho.  1 ,  pp.  17-25,  1979. 

2.  K.K.  Kuo,  J.H.  Koo,  T.R .  Davis,  G.R.  Coates,  " Transient  Combustion 
in  Mobile  Gas-Permeable  Propellants ,"  Acta  Astronautica,  Vol.  3, 
pp.  573-591 ,  1976. 

3.  E.B.  Fisher,  K.W.  Graves,  A.P.  Trippe,  " Application  of  a  Flame 
Spread  Model  to  Design  Problems  in  the  255  mm  Propelling  Charge ," 
12th  JAHHAF  Combustion  Meeting,  CPIA  Pub.  273,  Vol.  I,  pp.  199-219, 
December  1975. 

4.  H.  Krier,  S.S.  Kokhale,  " Modeling  of  Convective  Mode  Combustion 
Through  Granulated  Propellant  to  Predict  Detonation  Transition ," 

AJAA  Journal,  Vol.  16,  Ho.  2,  pp.  177-183 ,  1978. 

5.  P.S.  Gough,  "A  Two-Dimensional  Model  of  the  Interior  Ballistics  of 
Bagged  Artillery  Charges, "  Paul  Gough  Associcates ,  Inc.,  Report 
PGA-TR-81-1 ,  Portsmouth,  HH ,  1981. 
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analysis  for  the  ignition  phase  of  the  ballistic  cycle.  The  present 
analysis  and  computer  code  ALPHA6’  were  developed  by  Scientific 
Research  Associates  under  Ballistic  Research  Laboratory  sponsorship. 
This  computer  code  is  the  first  two-phase,  two-dimensional,  Navier- 
Stokes  analysis  to  simulate  the  interior  ballistic  cycle.  A  previous 
paper8  presented  the  initial  one-phase  flow  applications  of  the  ALPHA 
code  to  the  ballistic  cycle.  These  included  laminar  flow-adiabatic 
wall,  turbulent  flow-adiabatic  wall  and  laminar  flow-isothermal  wall 
simulations . 

This  report  presents  the  two-phase,  turbulent  flow  conservation 
equations  (References  6  and  7)  and  the  initial  two-dimensional 
two-phase,  nonreacting  viscous  laminar  flow  simulation  for  an  idealized 
ballistic  cycle.  These  two-phase  results  are  compared  with  a  one-phase 
simulation  of  the  cycle  to  illustrate  the  influence  of  solid  particles 
on  the  flow  development  in  the  tube. 


II.  GOVERNING  EQUATIONS 

The  governing  equations  for  a  two-phase,  two-dimensional  flow  in  a 
gun  tube  have  been  presented  previously  in  References  6  and  7.  Since 
the  large  number  of  propellant  particles  present  precludes  treating  each 
individually,  the  present  formulation  considers  averaged  flow  prop¬ 
erties.  These  equations  (References  6  and  7)  were  obtained  using  the 
formal  averaging  procedure  presented  by  Gough9’10  and 
Gough  and  Zwarts  (Reference  1),  with  extensive  reference  to  the  time- 

S~.  TF.  J~.  G%SeVCng~,  R.C.  Buggeln  and  H.  McDonald ,  "Development  of  a 
Two-Dimensional  Implicit  Interior  Ballistics  Code,"  US  Army 
ARRADCOM/Ballistic  Research  Laboratory,  Contract  Report 
ARBRL-CR-00411 ,  Aberdeen  Proving  Ground,  MD ,  January  1980  (AD  A084092). 

7.  H.J.  Gibeling  and  H.  McDonald,  " Development  of  a  Two-Dimensional 
Implicit  Interior  Ballistics  Code,"  US  Army  ARRADCOM/Ballistic 
Research  Laboratory,  Contractor  Report  ARBRL-CR-004S1 ,  Aberdeen 
Proving  Ground,  MD,  March  1981  (AD  A100276). 

8.  J.A.  Schmitt,  N.E.  Banks,  C.K.  Zoltani,  H.J.  Gibeling  end  T.L.  Mann, 
"Two-Phase  Viscous  Flow  Modeling  of  Interior  Ballistics ,  Algorithm 
and  Numerical  Predictions  for  an  Idealized  Lagrange  Gun," 

Proceedings  of  the  ASMS  Symposium  on  Computers  in  Flow  Predictions 
and  Fluid  Dynamic  Experiments ,  Edited  by  K.N.  Ghia,  T.J.  Mueller  and 
B.R.  Patel,  pp.  181-190,  1981. 

9.  P.S.,  Gough ,  " Derivation  of  Balance  Equations  for  Heterogeneous 
Two-Phase  Flow  by  Formal  Averaging, "  ARO  Workshop  on  Multiphase 
Flows,  USA  Ballistic  Research  Laboratory ,  pp.  71-80,  February  1978. 

10.  P.S.  Gough ,  "The  Flow  of  a  Compressible  Gas  Through  an  Aggregate  of 
Mobile,  Reacting  Particles ,"  Ph.D  Thesis,  Department  of  Mechanical 
Engineering ,  McGill  University ,  Montreal ,  1974. 
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averaging  approach  of  Ishii.11  In  this  approach  an  averaged  variable  is 
the  space-time  integral  of  the  product  of  the  local  variable  and  a 
weighting  function  which  reflects  the  influence  of  remote  points  on  the 
average  values.  The  local  equations  are  replaced  by  a  set  of  averaged 
equations  in  the  averaged  variables  for  each  spatial  position  and  time. 
The  two  phases  are  treated  distinctly  in  this  so-called  two-fluid  model 
with  momentum,  mass  and  energy  exchange  terms  between  the  two  phases 
being  replaced  by  appropriate  constitutive  relations.  The  governing 
partial  differential  equations  and  constitutive  relations  used  in  the 
present  analysis  are  given  in  the  Appendix. 

The  governing  equations  in  the  ALPHA  code  were  formulated  in 
conservation  form  by  application  of  a  Jacobian  transformation  to  the 
equations  in  cylindrical— polar  coordinates  (Reference  7).  A  consistent 
technique  for  determination  of  the  local  time-dependent  Jacobian  deter¬ 
minant  of  the  coordinate  transformation  was  implemented  to  reduce  geo¬ 
metrical  errors  in  computed  results.  An  axisymmetric  time-dependent 
adaptive  coordinate  system  for  interior  ballistics  flow  field  calcu¬ 
lations  is  utilized,  with  the  projectile  and  any  distinct  filler  ele¬ 
ments  treated  using  a  quasi-one-dimensional  lumped  parameter  analysis. 
The  governing  equations,  constitutive  relations,  and  the  time-dependent 
coordinate  system  have  been  incorporated  into  an  existing  computer  code 
which  solves  the  multidimensional  time-dependent  compressible  Navier- 
Stokes  equations  using  a  consistently  split,  linearized j ?bl^ck-implicit 
(LBI)  numerical  scheme  developed  by  Briley  and  McDonald.  " ’  Since 

the  LBI  algorithm  is  well  documented  in  References  12  and  13,  it  will 
not  be  discussed  further  here. 


III.  TWO-PHASE  FLOW  MODEL  PROBLEM 


Because  of  the  complex  nature  of  the  multidimensional,  two-phase 
flow  phenomena  occurring  in  a  gun  tube,  an  idealized  two-phase  flow  in  a 
Lagrange  gun  has  been  simulated.  This  problem  permits  verification  of 
the  two-phase  equation  system  in  the  ALPHA  code  and  provides  an 
understanding  of  some  basic  two-phase  flow  mechanisms  which  occur  in  a 
real  gun.  The  Lagrange  gun  is  a  smooth  tube  of  constant  radius  which  is 
closed  at  one  end  by  the  breech.  The  model  combustion  chamber  is 
bounded  by  the  breech,  a  flat  based  projectile,  and  the  tube  wall.  The 
chamber  is  filled  initially  with  a  high  pressure,  high  temperature 
gas-particle  mixture  with  a  small  gas-phase  region  adjacent  to  the 


11.  M.  Ishii,  Theimio-Fluid  Dynamic  Theory  of  Tuo^Phase  Plow., 

Eyrolles,  Paris,  1975. 

12.  W.E.  Briley  and  H.  McDonald,  "Solution  of  the  Multidimensional 
Compressible  Navier-Stokes  Equations  by  a  Generalized  Implicit 
Method,"  Journal  of  Computational  Physics,  Vol.  24,  No.  4,  pp. 
372-397,  1977. 

13.  W.R.  Briley  and  H.  McDonald ,  "On  the  Structure  and  Use  of  Linearized 
Block  Implicit  Schemes,"  Journal  of  Computational  Physics^,  Vol.  34, 
No.  1,  pp.  54-73,  1980. 
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projectile  base.  In  the  present  case,  the  particles  were  chosen  to  be 
spherical  with  a  constant  0.1  ram  diameter  throughout.  No  combustion  of 
the  particles  is  permitted  in  this  calculation,  so  that  this  simulation 
of  the  ballistic  cycle  approximately  represents  that  of  a  real  gun  if 
the  propellant  combustion  ceases  before  the  projectile  begins  moving  at 
time  t  *  0.  As  in  prior  one-phase  Lagrange  gun  simulations  (Reference 
8),  frictional  forces  between  the  projectile  and  the  tube  wall  are  ne¬ 
glected.  For  this  initial  two-phase  simulation  the  flow  is  assumed  to 
remain  laminar  throughout  the  cycle,  and  an  adiabatic  wall  condition  is 
applied.  Even  though  the  actual  flow  in  a  gun  tube  is  turbulent  and  the 
heat  transfer  to  the  walls  influences  the  flow  field,  the  present  simu¬ 
lation  provides  insight  into  the  effects  of  solid  particles  on  the  vis¬ 
cous  flow  development  in  the  tube.  The  required  governing  equations  are 
the  gas  and  solid  phase  continuity  and  momentum  conservation  equations, 
gas  phase  energy  conservation  equation  and  the  solid  phase  heat  conduc¬ 
tion  equation.  Heat  transfer  between  the  gas  and  solid  phases  is  accom¬ 
plished  by  the  interphase  heat  transfer  coefficient  model  with  the  par¬ 
ticle  surface  temperature  advanced  in  time  after  the  coupled  solution  of 
the  gas  and  solid  phase  continuity  and  momentum  equations  and  the  gas 
phase  energy  equation.  However,  the  test  case  was  run  with  thermally 
insulated  particles,  i.e.,  no  heat  transfer  between  the  gas  and  solid 
phase  was  permitted  and  the  particle  internal  energy  remained  constant 
throughout  the  cycle.  The  particle  size  is  assumed  to  remain  constant, 
and  the  gas  phase  specific  heat  (at  constant  pressure)  and  mixture  mo¬ 
lecular  weight  are  also  assumed  to  remain  constant  throughout  the  cycle. 
The  Nobel -Abel  equation  of  state,  Sutherland's  laws  for  viscosity  and 
thermal  conductivity,  and  the  Fourier  heat  conduction  relation  are 
employed  for  the  gas  phase. 

The  following  coupled  boundary  conditions  are  utilized  for  the 
two-phase  calculation:  at  the  breech,  projectile  and  tube  wall,  no-slip 
conditions  are  imposed  on  the  gas  phase  velocities  and  the  gas  phase 
density  is  determined  from  the  normal  momentum  equation.  At  the  solid 
surfaces,  the  adiabatic  condition  of  zero  normal  derivative  of  gas  tem¬ 
perature  is  applied.  At  the  breech  and  tube  wall,  the  normal  component 
of  solid  phase  velocity  is  assumed  to  be  zero,  while  the  tangential  com¬ 
ponent  is  determined  by  assuming  that  there  is  no  solid  phase  momentum 
decrease  due  to  interaction  with  the  wall.  This  latter  condition  is  ap¬ 
proximated  by  assuming  a  zero  normal  derivative  of  the  solid  phase  tan¬ 
gential  velocity  component.  The  solid  phase  partial  density  (related  to 
void  fraction)  at  the  breech  is  obtained  by  application  of  the  solid 
phase  continuity  equation  and  at  the  tube  wall  is  obtained  by  assuming  a 
zero  normal  derivative  of  solid  phase  partial  density.  At  the  projec¬ 
tile  base,  the  solid  phase  velocity  and  partial  density  are  set  to  zero 
(unit  void  fraction)  since  there  is  only  gas  present  at  that  location. 

At  the  tube  centerline,  the  symmetry  conditions  of  zero  radial  veloc¬ 
ities  and  zero  radial  derivatives  of  all  other  variables  are  imposed. 

The  computational  mesh  employed  In  the  calculations  consists  of  61 
uniformly  spaced  mesh  points  in  the  axial  direction  and  19  nonuniformly 
spaced  mesh  points  in  the  radial  direction.  The  mesh  points  are  con¬ 
centrated  near  the  tube  wall  in  the  radial  direction  In  order  to 
resolve  the  gas  phase  boundary  layer;  the  first  point  is  at  the  tube 
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centerline  (10  ms  from  the  wall),  the  tenth  is  0.306  mm  from  the  wall, 
the  seventeenth  is  7.56  ym  from  the  wall,  and  the  nineteenth  point  is 
at  the  wall.  The  gun  tube  configuration,  gas  and  solid  phase 
properties,  and  initial  conditions  for  the  calculation  are  given  in 
Table  1. 

TABLE  1.  TWO-PHASE  LAGRANGE  GUN  PARAMETERS 


Bore  Diameter 

20mm 

Combustion  Chamber  Length 

0.175m 

Total  Bore  Length 

1.290m 

Projectile  Mass 

0.120kg 

Gas  Molar  Mass 

23. 8g/mole 

Covolume 

1.08  •  10" 

Ratio  of  Cas  Specific  Heats 

1.271 

Solid  Particle  Diameter 

0.1mm 

Particle  Density 

1500  kg/m3 
£ 

Initial  Gas  Pressure 

300  •  10  1 

Initial  Gas  Temperature 

3000  K 

Initial  Void  Fraction;  0  -  Z  -  0.169m 

0.94 

0.169  <  Z  <  0.175m 

1.0 
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/kg 


A  one-phase  flow  Lagrange  gun  simulation  was  computed  with  nearly 
the  same  computational  mesh  ( 19  x  49)  for  comparison  with  the  two-phase 
calculation.  The  additional  grid  points  in  the  axial  direction  were 
used  in  the  two-phase  calculation  because  of  the  large  void  fraction 
gradient  occurring  near  the  projectile  base.  The  one-phase  simulation 
has  the  same  parameters  as  Table  I  except,  of  course,  that  the  void 
fraction  is  unity  throughout.  As  previously  reported  by  Schmitt  et  al 
(Reference  8),  the  one-phase  Lagrange  gun  results  compare  favorably  with 
both  a  one-dimensional  analytic  isentropic  solution1  which  is  valid  for 
early  times  in  the  cycle  and  the  numerical  calculation  of  Heiser  and 
Hensel.1^  A  one-phase  calculation  run  with  the  present  computer  code 
was  compared  to  an  identical  calculation  run  with  the  code  utilized  by 
Reference  8  in  order  to  verify  the  present  one-phase  results.  These 
results  give  confidence  that  the  numerical  procedure  is  capable  of 
capturing  the  pressure  waves  propagating  in  the  gun  tube  as  long  as  the 
spatial  and  temporal  resolution  remain  adequate. 

The  two-phase  calculation  was  started  with  a  time  step  of  2.0  ps, 
and  the  time  step  was  permitted  to  increase  to  maintain  a  predetermined 
maximum  change  in  the  solution  over  each  time  step.  The  maximum  time 
step  reached  in  this  case  was  about  4.3  ps.  The  one-phase  calculation 
was  run  with  a  constant  time  step  of  5.0  ps,  so  that  temporal  accuracy 
in  the  two  cases  should  be  comparable.  The  calculations  were  both  ter¬ 
minated  at  muzzle  exit  time.  The  results  are  compared  shortly  before 
muzzle  exit  time  with  the  projectile  speed,  position,  and  problem  time 
of  581  m/s,  1.282  m,  and  2.575  ms,  respectively,  for  the  one-phase  simu¬ 
lation  and  559  m/s,  1.289  m,  and  2.673  ms,  respectively,  for  two-phase 
simulation.  The  one-phase  calculation  was  actually  run  in  the  two-phase 
mode  with  near  unity  void  fraction  (0.999999),  but  solving  only  the 
gas-phase  conservation  equations.  The  one-phase  calculation  required 
518  time  steps  and  0.82  hours  of  CP  time  on  the  BRL  CDC  7600  with  the  19 
x  49  mesh  and  4  coupled  equations;  the  two-phase  calculation  required 
891  time  steps  and  3.06  hours  of  CP  time  with  the  19  x  61  mesh  and  7 
coupled  equations • 

The  comparison  of  the  boundary  layer  development  in  the  one-phase 
and  two-phase  cases  shows  some  interesting  results.  The  boundary  layer 
axial  velocity  profiles  are  shown  in  Figures  1  through  4,  for  several 
axial  positions,  from  which  it  is  seen  that  the  gas  velocities  are  smal¬ 
ler  in  the  two-phase  case  than  in  the  one-phase  case,  reflecting  the 
energy  expended  to  accelerate  the  solid  particles.  The  solid  phase  lags 
behind  the  gas  in  the  core  flow  as  expected;  the  particle  velocity  de¬ 
creases  very  little  near  the  wall  because  of  the  relatively  large  particle 

14.  E.H.  Love  and  F.B.  Pidduck,  "Lagrange9 8  Ballistic  Problem ,"  Phil . 
Trane.  Roy.  Soc.,  Vol.  222 ,  pp.  167-226,  1921-22. 

15.  R.  Reiser  and  D.  Hensel ,  " Berechnung  der  Gastromung  in  einem 
Waffenrohr  mit  Hilfe  des  Zweidimensionalen  AMI-Modells  ,"  El/81, 
Emest-Mach-Institut ,  Abteilung  fur  Ballistik,  Weil  am  Rhein,  West 
Germany,  January  1981. 
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Distance  from  wall  (mm) 


gas,  one-phase  calculation,  t  =  2.575  ras 
gas*  two-phase  calculation,  t  *  2.673  ras 
particle,  two-phase  calculation,  t  =  2.673  ms 


0  200  400  600 


Axial  Velocity  (m/s) 

Figure  1  -  Comparison  of  one-  and  two-phase  tube  wall  boundary 
layer  profiles  at  muzzle  exit  time  for  Z  =  0.43  ra. 
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Distance  from  wall  (mm) 


gas,  one-phase  calculation,  t  =  2.575  ms 
gas,  two-phase  calculation,  t  =  2.673  ms 
particle,  two-phase  calculation,  t  =  2.673  ms 


Axial  Velocity  (m/s) 

Figure  2  -  Comparison  of  one-  and  two-phase  tube  wall  boundary 
layer  profiles  at  muzzle  exit  time  for  Z  s  0,64  in. 
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Distance  from  wall  (mm) 


gas,  one-phase  calculation,  t  =  2.575  ms 
gas,  two-phase  calculation,  t  =  2.673  ms 
particle,  two-phase  calculation,  t  -  2.673  ms 


Axial  Velocity  (m/s) 

Figure  3  -  Comparison  of  one-  and  two-phase  tube  wall  boundary 
layer  profiles  at  muzzle  exit  time  for  Z  -  0.86  m. 
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Distance  from  wall  (mm) 


gas,  one-phase  calculation,  t  =  2.575  ms 
gas,  two-phase  calculation,  t  =  2.673  ms 
particle,  two-phase  calculation,  t  =  2.673  ms 


0.64 


0.32 


0 

0  200  400  600 

Axial  Velocity  (m/s) 

Figure  4  -  Comparison  of  one-  and  two-phase  tube  wall  boundary 
layer  profiles  at  muzzle  exit  time  for  Z  =  1.07  ra. 
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size  and  the  assumption  of  no  particle  momentum  loss  at  the  wall. 

The  transverse  velocity  profiles  across  the  entire  tube  are  shown  in 
Figures  5  through  8  for  the  same  axial  positions  as  in  Figures  1  through 
4.  The  transverse  velocities  are  quite  small  compared  to  the  axial  vel¬ 
ocities  for  this  relatively  simple  case.  Also,  the  velocities  near  the 
juncture  of  the  tube  wall  and  the  projectile  base  are  strongly  influ¬ 
enced  by  the  treatment  of  that  point.  In  the  present  calculations  a 
perfect  seal  was  assumed  so  that  this  corner  juncture  is  a  singular 
point.  In  order  to  obtain  more  accurate  predictions  in  the  corner 
region,  the  gap  region  between  the  tube  wall  and  the  projectile  may  have 
to  be  included  in  the  computational  domain.  The  positive  radial  veloci¬ 
ties  (toward  the  tube  wall)  observed  near  the  projectile  base  (Figure  8) 
are  expected  on  physical  grounds.  However,  it  has  been  observed  that 
the  quantitative  results  near  the  projectile  base  are  sensitive  to  the 
axial  grid  spacing,16  so  that  a  locally  refined  mesh  should  be  utilized 
in  that  region. 


The  character  of  the  tube  wall  boundary-layer  is  noticeably  dif¬ 
ferent  in  the  one-  and  two-phase  cases,  as  seen  in  Figures  1  through  8. 
The  axial  distribution  of  boundary  layer  displacement  thickness  (Figure 
9)  gives  another  indication  of  the  influence  of  solid  particles  on  the 
wall  boundary  layer.  Here,  the  displacement  thickness  (6*)  is 


s*(z,n  = 


yotr ,  z,  1 )  w  ( r ,  2,  1)  1  r 

pc(  2,1)  wc(^  iT~  J  *R  d( 


(3) 


where  Che  subscript  c  denotes  the  centerline  value,  w  is  the  axial  gas 
velocity,  p  is  the  gas  density,  and  R  is  the  tube  radius.  It  can  be 
seen  that  the  displacement  thickness  is  larger  in  the  two-phase  case  for 
most  of  the  distance  from  the  projectile  to  the  breech.  Furthermore, 
the  wall  shear  stress  is  larger  in  the  two-phase  case.  If  the  thermal 
boundary  layer  for  a  flow  with  wall  heat  transfer  exhibits  similar  be¬ 
havior,  an  increase  in  wall  heat  flux  would  be  expected  due  to  the 
presence  of  solid  particles  of  the  present  size. 


The  void  fraction  distribution  shown  in  Figure  10  illustrates  the 
increase  in  thickness  of  the  gas  zone  (a  *  1.0)  near  the  projectile; 
however,  this  is  essentially  a  one-dimensional  process  since  no  particle 
momentum  loss  at  the  wall  has  been  assumed  and  the  particle  size  is  the 
same  order  as  the  boundary  layer  thickness.  Since  the  particle  size  is 
comparable  to  the  gas  boundary  layer  thickness,  the  interphase  correla¬ 
tions  must  be  examined  to  include  near-wall  effects  and  a  particle-wall 
interaction  model  should  be  developed  to  account  for  particle  momentum 
transfer  to  the  wall.  The  axial  distribution  of  velocities  for  the  one- 
and  two-phase  calculations  are  compared  in  Figures  11  through  13.  The 

16.  J.A .  Schmitt  and  T.L.  Mann>  "Calculation  of  the  Compressible  Flow  in 
the  Lagrange  Gun  by  the  Interior  Ballistics  Algorithm  Alpha /'  Pro¬ 
ceedings  of  the  DEA-G-1060  Meeting 3  Eglin  Air  Force  Base3  October 
1980. 
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Distance  from  wall  (mm) 


gas,  one-phase  calculation,  t  =  2.575  ms 

gas,  two-phase  calculation,  t  =  2.673  ms 
particle,  two-phase  calculation,  t  =  2.673 


Radial  Velocity  (m/s) 

Figure  5  -  Comparison  of  one-  and  two-phase  radial  velocity 
profiles  at  muzzle  exit  time  for  Z  =  0.43  m. 
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ms 


Distance  from  wall  (mm) 


— 


gas,  one-phase  calculation,  t  *  2.575  ms 
gas,  two-  phase  calculation,  t  =  2.673  ms 
particle,  two-phase  calculation,  t  =  2.673  ms 


Figure  6  -  Comparison  of  one-  and  two-phase  radial  velocity 
profiles  at  muzzle  exit  time  for  Z  =  0.64  m. 
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Distance  from  wall  (mm) 


gas,  one-phase  calculation,  t  =  2.575  ms 
gas,  two-phase  calculation,  t  *  2.673  ms 
particle,  two-phase  calculation,  t  =  2.673  ms 


Figure  7  -  Comparison  of  one-  and  two-phase  radial  velocity 
profiles  at  muzzle  exit  time  for  Z  =  0.86  m. 
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Distance  from  wall  (mm) 


gas,  one-phase  calculation,  t  =  2.575  ms 
gas,  two-phase  calculation,  t  =  2.673  ms 
particle,  two-phase  calculation,  t  =>  2.673  ms 


Figure  8  -  Comparison  of  one-  and  two-phase  radial  velocity 
profiles  at  muzzle  exit  time  for  Z  =  1.07  m. 
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Displacement  thickness  (mm) 


one-phase  calculation,  t  =  2.575  ms 
two-phase  calculation,  t  =  2.673  ms 


Distance  from  breech  (m) 

Figure  9  -  Comparison  of  boundary  layer  displacement  thickness 

for  one-  and  two-phase  calculations  at  muzzle  exit  time. 
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Distance  from  breech  (m) 

Figure  10  -  Axial  distribution  of  void  fraction  at  tube 
centerline  at  muzzle  exit  time,  t  =  2.673  ms. 
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Axial  velocity  (m/s) 


- -  gas,  one-phase  calculation,  t  =  2.575  ms 

-  8as>  two-phase  calculation,  t  =  2.673  ms 


particle,  two-phase  calculation,  t  =  2.673  ms 


Distance  from  breech  (m) 

Figure  11  —  Comparison  of  axial  velocity  distributions  at  tube 
centerline  for  one-  and  two-phase  calculations  at 
muzzle  exit  time. 


24 


Axial  velocity  (m/s) 


gas,  one-phase  calculation,  t 
gas,  two-phase  calculation,  t 


2.550  ms 
2.635  ms 


particle,  two-phase  calculation,  t  =  2.635  ms 


Distance  from  breech  (m) 

Figure  12  -  Comparison  of  axial  velocity  distributions 
at  r  =  9.9924  mm  for  one-  and  two-phase 
calculations  near  muzzle  exit  time. 
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Radial  velocity  (m/s) 


gas, 

gas. 


one-phase  calculation, 
two-phase  calculation. 


t 

t 


2 . 550  ms 
2.635  ms 


particle,  two-phase  calculation,  t  =  2.635  ms 


Distance  from  breech  (m) 

Figure  13  -  Comparison  of  radial  velocity  distributions  at 

r  =  9.9924  mm  for  one-  and  two-phase  calculations 
near  muzzle  exit  time. 


lower  gas  velocity  in  the  two-phase  results  reflects  the  energy  required 
to  accelerate  the  solid  particles.  Radial  direction  distributions  of 
pressure,  density  and  temperature  are  shown  in  Figures  14,  15  and  16, 
respectively,  while  axial  direction  distributions  are  shown  in  Figures 
17,  18  and  19,  respectively.  An  abrupt  change  in  pressure  near  the  tube 
wall-projectile  corner  (Figure  17)  is  evident  due  to  the  discontinuous 
change  in  velocity  from  the  projectile  base  to  the  tube  wall.  The 
density  and  temperature  variation  near  the  projectile  base  are  known  to 
be  inaccurate  with  a  uniform  axial  mesh  (Reference  16),  One-phase 
calculations  performed  by  Reference  16  with  a  refined  axial  mesh  near 
the  projectile  base  did  not  exhibit  either  the  increase  in  temperature 
or  decrease  in  density  as  the  projectile  is  approached  near  the  tube 
wall.  As  mentioned  previously,  this  problem  area  must  be  addressed 
carefully  to  improve  computed  results  in  the  corner  region. 

In  the  present  calculations,  total  system  mass  and  energy  are 
computed  periodically  as  a  check  on  the  spatial  and  temporal  accuracy  of 
the  finite  difference  computations.  In  the  one-phase  case,  the  total 
energy  of  the  system  consists  of  gas  internal  energy,  gas  kinetic  energy 
and  projectile  kinetic  energy.  For  the  one-phase  results,  the  final 
mass  and  energy  of  the  system  differed  from  starting  values  by  -0.16% 
and  -0.23%,  respectively.  In  the  two-phase  case,  total  energy  consists 
of  gas  internal  energy,  gas  kinetic  energy,  particle  internal  energy 
(constant  for  thermally  insulated  particles),  particle  kinetic  energy, 
and  projectile  kinetic  energy.  For  the  two-phase  results,  the  final 
mass  and  energy  of  the  system  differed  from  starting  values  by  0.17%  and 
-3.11%,  respectively.  The  larger  discrepancy  in  energy  conservation  in 
the  two-phase  calculation  is  believed  to  be  associated  with  inadequate 
resolution  of  the  interface  between  the  pure  gas  zone  and  the 
gas-particle  mixture  near  the  projectile.  Ultimately,  an  adaptive  mesh 
transformation  should  be  implemented  to  concentrate  more  grid  points 
near  such  an  interface  in  order  to  reduce  local  spatial  truncation 
errors . 


IV.  CONCLUSIONS 

The  present  two-phase,  averaged  Navier-Stokes  analysis  has  been 
applied  to  the  two-phase  flow  in  a  Lagrange  gun.  A  consistently  split, 
linearized,  block  implicit  numerical  procedure  has  been  employed  to 
solve  the  governing  partial  differential  equations.  This  initial 
two-phase  flow  application  has  been  restricted  to  laminar  flow  and 
adiabatic  wall  conditions;  however,  the  ability  of  the  computer  code  to 
solve  the  governing  equations  has  been  demonstrated. 

It  has  been  noted  that  the  interphase  correlations  must  be  examined 
to  include  near  wall  effects,  and  the  influence  of  particle-wall 
interaction  must  be  determined.  These  areas  will  require  more  detailed 
experimental  measurements  than  are  presently  available,  and  to  this  end 
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Distance  from  wall  (mm) 


-  gas,  one-phase  calculation,  t  =  2.550  ms 

-  —  —  gas,  two-phase  calculation,  t  =  2.635  ms 


Pressure  (10^  Pa) 


Figure  14  -  Comparison  of  boundary  layer  pressure  profiles 

for  one-  and  two-phase  calculations  at  two  .axial 
locations  near  muzzle  exit  time. 
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-  —  gas,  two-phase  calculation,  t 


2.550  ms 
2.635  ms 


23.0  28.0  33.0 


3 

Gas  density  (kg/m  ) 

Figure  15  -  Comparison  of  boundary  layer  density  profiles 
for  one-  and  two-phase  calculations  at  two 
axial  locations  near  muzzle  exit  time. 
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gas,  one-phase  calculation,  t 
gas,  two-phase  calculation,  t 


2.550  ms 
2.635  ms 


Gas  temperature  (K) 

Figure  16  -  Comparison  of  boundary  layer  temperature 

profiles  for  one-  and  two-phase  calculations 
at  two  axial  locations  near  muzzle  exit  time. 
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Figure  17 
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gas,  two-phase  calculation,  t 
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-  Comparison  of  pressure  distributions  for  one-  and 
two-phase  calculations  at  two  radial  locations 
near  muzzle  exit  time. 
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Gas  density  (kg/m  ) 


gas,  one-phase  calculation,  t  =  2.550  ms 
gas,  two-phase  calculation,  t  =  2.635  ms 


Distance  from  breech  (m) 

Figure  18  -  Comparison  of  density  distributions  for  one-  and 
two-phase  calculations  at  two  radial  locations 
near  muzzle  exit  time. 
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-  gas,  one-phase  calculation,  t  =  2.550  ms 
Sas»  two-phase  calculation,  t  =  2.635  ms 


Distance  from  breech  (m) 

Figure  19  -  Comparison  of  temperature  distributions  for  one- 

and  two-phase  calculations  at  two  radial  locations 
near  muzzle  exit  time. 
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benchmark  experiments  in  a  simulated  ballistic  environment  have  been 
designed  by  BRL  personnel  and  J.  Whitelaw  of  Imperial  College.  The 
resulting  data  will  be  used  in  future  modeling  efforts  for  both 
improvement  of  the  interphase  correlations  and  verification  of 
turbulence  models  for  two-phase  flows  in  gun  tubes. 
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APPENDIX  A 


GOVERNING  EQUATIONS 


A.  Governing  Equations 


The  averaged  time-dependent  governing  partial  differential  equa¬ 
tions  for  the  turbulent,  two-phase  flow  in  a  gun  tube  were  presented  in 
References  6  nnd  7.  The  equations  are  summarized  here  in  dimen¬ 
sional  form  with  corrections  where  required. 

The  gas  and  solid  phase  continuity  equations  are,  respectively. 


d(ap ) 

ai“ 


V  •  (  apu)  =  F 


(A-l ) 


d  ( I  -  a  )pp 
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where  the  void  fraction  (or  porosity)  a  is  the  ratio  of  volume  occupied 
by  the  gas  phase  to  the  total  volume.  The  averaged  densities  and  vel¬ 
ocities  of  the  gas  and  solid  phases  are  p,  U,  pp,  Up,  respectively. 

The  interphase  mass  transfer  rate  T ,  Eq .  (A-12),  is  due  to  propel¬ 
lant  burning.  The  gas  and  solid  phase  momentum  equations  are,  respec¬ 
tively, 
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where  p  is  the  averaged  pressure.  The  tensors  it  and  tt^,  Equations 
(A-16)  and  (A-19),  are  the  averaged  molecular  and  turbulent  stress  ten¬ 
sors  in  the  gas  phase*  respectively.  The  constitutive  relations  for  the 
propellant  grain  surface  area  Sp,  volume  Vp,  the  interphase  drag 
<F> ,  and  the  isotropic  intergranular  stress  Rp  are  given  by  Equations 
( A-13) ,  (A-14),  (A-22) ,  and  (A-25),  respectively. 

The  gas  phase  energy  equation  is 


d(aph)  _  r  _ 

“ 0 ,  —  4  V-fa^pDh)  =  -  V  •  [  a  (  q  4  q 
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where  h  is  the  averaged  static  enthalpy  of  the  gas.  The  mean  flow 
energy  dissipation  $ ,  the  turbulent  energy  dissipation  t^ie  inter¬ 

phase  energy  transfer  A,  and  the  molecular  and  turbulent  heat  flux  vec¬ 
tors  2j  and  are  given  by  Equations  (A-26),  (A-27),  (A-28) ,  (A-34) 
and  (A-35),  respectively.  Also,  the  symbol  D/Dt  denotes  the  substantial 
derivative,  i.e., 


D 

Dt 


71  +  uv 


(A-6) 


The  Nobel -Abel  equation  of  state,  Equation  (A-37),  relates  the  gas  pres¬ 
sure  to  the  density  and  temperature.  The  turbulence  kinetic  energy 
equation  is  given  by 


d(a  />k)  u. 
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dt 
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4  S. 


where  =  1.0  is  used,  and  the  interphase  turbulence  production, 
Sfc,  is  neglected  at  present.  The  rate  of  strain  tensor  D  and  the 
turbulent  viscosity  pT  are  given  by  Equations  (A-17)  and  (A-20) , 
respectively. 
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Exclusion  of  gas  phase  chemical  reactions  from  the  present  analysis 
eliminates  the  need  for  gas  phase  and  gasified  propellant  species  mass 
fraction  equations.  However,  in  order  to  consider  several  types  of  pro¬ 
pellants  in  the  system,  transport  equations  for  the  reciprocal  gas  mix¬ 
ture  molecular  weight  parameter  Z  =  Ry/Wgj  and  the  specific  heat  at 
constant  pressure  must  be  solved;  these  are 


and 


d\ 


(apZ)  +  V-  [ap  U Z)  =  V-  [aTmVz]  +  2  T  (A-8) 
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where  Zp  and  (cp)p  are  the  reciprocal  molecular  weight  parameter 
and  specific  heat  of  the  propellant,  respectively,  and  rm  * 

Peff/Sceff  =  heff  =  ^  +  ^T,  and  the  effective  Schmidt  number  is 
taken  as  Sceff  =  0.9. 

At  present,  the  propellant  grains  are  assumed  to  be  spherical.  The 
surface  temperature  of  these  particles  is  desired  to  determine  ignition, 
the  rate  of  heat  transfer  between  the  gas  and  solid  phases,  and  possibly 
the  propellant  burning  rate.  By  assuming  that  the  penetration  depth  of 
the  thermal  wave  into  the  propellant  grains  is  small  compared  to  the 
grain  dimensions,  it  is  permissible  to  use  a  one-dimensional  approxi¬ 
mation  to  obtain  the  propellant  surface  temperature.  Following  the 
motion  of  a  given  particle,  the  heat  conduction  equation  for  the  propel¬ 
lant  temperature  Tp  is 


(A-10) 


where  Tp  =  Tp(r;x,t)  is  the  phase-averaged  temperature  within  a 
representative  particle,  r  is  the  radial  location  within  the  spherical 
particle,  dp  is  the  thermal  diffusivity  of  the  particles,  and 
(d/dt)r  denotes  the  Lagrangian  time  derivative  at  constant  r  within 
the  particle.  The  relevant  boundary  condition  for  Equation  (A-10)  is  a 
specified  surface  heat  flux.  Equation  (A- 41),  which  is  determined  using 
the  interphase  heat  transfer  relation,  Equation  (A—  29) ,  and  a  heat  feed¬ 
back  relation,  Equation  (A- 42),  identified  by  Gough. A-1  The  actual 

A-l.  P.S.  Gough ,  "Numerical  Analysis  of  a  Two-Phase  Flow  with  Explicit 
Internal  Boundaries,"  IHCFt  77- 5,  Naval  Ordnance  Station,  Indian 
Head,  MD ,  April  1977. 
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burning  is  modeled  as  a  regression  of  the  propellant  surface  with  a 
resulting  deposition  of  mass,  momentum  and  energy  into  the  gas  phase, 
The  equation  for  average  particle  radius,  rp,  including  turbulent 
diffusion  is  r  ■* 


a["-a), Ppfpl  4?.[(l-a)PpUprp] 

•  v[(l-a)rmVrp]-(l-a)Pp(l  +  rp^) 


<  d> 


(A-ll) 


where  the  surface  regression  rate  <d>  is  given  by  Equation  (A-15).  The 
relation  for  T,  Equation  (A-12),  has  been  utilized  in  order  to  cast  the 
particle  radius  equation  into  the  above  form. 


The  influence  of  the  hot  primer  gas  may  be  included  in  the  pro¬ 
cedure  in  two  ways.  The  preferred  approach  is  the  application  of  mass 
flux  and  temperature  boundary  conditions  near  the  tube  centerline  for  a 
center  core  primer  or  at  the  breech  end  for  a  base  primer.  An  alterna¬ 
tive  method  would  be  the  introduction  of  mass,  momentum,  and  energy 
sources  into  the  governing  gas  phase  equations  to  represent  the  primer 
discharge . 

The  governing  Equations  (A-l)  through  (A-  9 )  and  (A-ll)  are  a 
system  of  general  time-dependent  nonlinear  partial  differential 
equations  for  two-phase  flow.  The  applicability  of  these  equations  to  a 
specific  two-phase  flow  depends  on  the  validity  of  both  the  averaging 
procedure  and  the  constitutive  relations,  Equations  (A-12)  through 
(A-40),  which  complete  the  analysis.  In  order  to  solve  these  equations 
they  must  be  written  in  a  specific  coordinate  system.  The  governing 
equations  in  the  ALPHA  code  were  formulated  in  conservation  form 
(Reference  7)  by  application  of  a  nonorthogonal  Jacobian  transformation 
to  the  Equations  (A-l)  through  (A-  9 )  and  (A-ll)  written  in  cylindrical 
polar  coordinates.  A  consistent  technique  for  determination  of  the 
local  time-dependent  Jacobian  determinant  of  the  coordinate  transfor¬ 
mation  was  implemented  in  order  to  reduce  geometrical  errors  in  the 
computations,  as  suggested  by  Thomas  and  Lombard. A- 2  resuiting 

equations  are  solved  for  the  following  dependent  variables:  the  radial, 
angular  and  axial  components  of  gas  velocity  ti,  the  gas  phase  partial 
density  ap ,  the  gas  phase  static  enthalpy  h,  specific  heat  Cp, 
reciprocal  molecular  weight  Z,  and  turbulence  kinetic  energy  k.  The 
solid  phase  dependent  variables  are  the  momentum  flux  components  of 
(l-a)ppUp,  the  solid  phase  partial  density  (l-a)pp,  and  the 
particle  radius  parameter  (l-a)pprp.  The  solid  particle  surface 
temperature  is  determined  by  solution  of  Equation  (A- 10)  using  boundary 
conditions  Equations  (A-41)  and  (A-42)  after  solution  of  the  fully 
coupled  system  Equations  (A-l)  through  (A-  9)  and  (A-ll).  Since  the 
solution  of  Equations  (A-l)  through  (A-  9)  and  (A-ll)  is  an  implicit 
one,  the  solid  particle  surface  temperature  must  be  treated  explicitly 

A-2.  P.D .  Thomas  and  C.K .  Lombard ,  "Geometric  Conservation  Law  and  Its 
Application  to  Flow  Computations  on  Moving  Grids,"  AIAA  Journal, 

Vol.  17,  No.  10,  p.  1030,  1979. 
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at  present.  This  may  result  in  a  time  step  limitation  to  preserve 
transient  accuracy;  however,  preliminary  results  indicate  that  this  is 
not  a  severe  restriction. 

B.  Constitutive  Relations 


Here  we  present  the  various  constitutive  relations  required  for 
closure  of  the  two-phase  flow  Equations  (A-l)  through  (A-ll).  The  gas 
phase  mass  source  due  to  propellant  burning  is 


-p  /  ,  \  P 

r  =  (  I  -  a )  -  —  <  d  > 

VP 


(A-12) 


where  the  grain  surface  area,  grain  volume,  and  the  steady  state  surface 
regression  rate  are  given  by 


vp  s  f  w,p 


(A-13) 


SP  -  rp 


<d>  =  Bf  +  B2pn 


(A-l A ) 
(A-l 5) 


respectively.  The  constants  Bj,B2  and  n  are  known  for  a  given 
propellant.  The  gas  molecular  stress  tensor  for  a  Newtonian  fluid  is 


7T  1  D  *  (-§-/*  "  Kb)  V-  U  I  (A'16> 

where  p  is  the  molecular  viscosity,  Kg  is  the  bulk,  viscosity,  and  I  is 
the  identity  tensor.  The  rate  of  strain  tensor  is 

—  I  r 

D  5  ~  [(VU)  4  (Vu)TRANSPOSEj  (A-17) 

and  the  molecular  viscosity  is  determined  from  Sutherland’s  law 


F  _  (  2  7o'isl  (a'18> 

Vo  '  V  T0  /  1+S, 
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where  Sj  =  110  K,  TQ  =  300  K,  and  PQ  =  1.8463* 10"5  Pa-s.  The 

bulk  viscosity  Kg  is  assumed  to  be  zero.  The  turbulent  stress  tensor 

is  modeled  using  an  isotropic  eddy  viscosity  formulation, 


ir7  =  -p  u'u'  =  2/i.jD  — —  (/iTV-U+/>k)I 


(A-19) 
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where  the  mixing  length,  £,  is  specified  algebraically  and  is  a 
function  of  the  turbulence  Reynolds  number,  RT=p^,/p. 


The  interphase  drag  relation  <F>  is  determined  by  a  combination  of 
Ergun's  relation^-^ f or  a  packed  bed  of  spheres  and  a  fit  to  the  drag 
coefficient  data  for  isolated  spheres.^-4 
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/1-3.  S~.  Ergun ,  rrFTuui  Flow  Through  Packed  Columns ,"  Chem.  Eng.  Progr.  , 

Vol.  48,  p.  89,  1952. 
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A.  Celmiy&,  Private  communication,  US  Army  Ballistic  Research 
Laboratory,  1982. 
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where  Ur  =  U  -  Up. 

The  intergranular  stress  relation  Rp  which  is  independent  of  the 
loading  history  is  given  by 
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where  ap  is  the  speed  of  sound  in  the  solid  phase  which  is  specified 
for  a  given  propellant,  and  ac  is  a  critical  porosity  above  which 
there  is  no  direct  contact  between  particles.  The  mean  flow  dissipation 
function  $  is  given  by 


^  "  2P  D  :  D  "  P  ~  Kb)(  V-U)? 


(A-26) 


and  the  turbulence  energy  dissipation  is  approximated  by 
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The  interfacial  energy  transfer  A  is  defined  by 
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where  hcomb  is  the  energy  released  per  unit  mass  due  to  combustion  of 
the  propellant.  The  interphase  heat  transfer  <q>  is  given  by 


<q>  »  h,(T-  Tps) 

where  the  total  heat  transfer  coefficient  ht 
h|  '  %  +fPtr(f  +T 


is 


(A-29) 


(A-30) 


where  Cp  is  the  propellant  emissivity,  o  is  the  Stefan-Boltzmann 
constant,  and  the  propellant  surface  temperature  is  determined  from  the 
solution  of  Equation  (A-10).  The  Nusselt  number  correlation  used  is 


Nup  =  2.0  +  0.4  Re*'3  Pf  1/3 


( A-31 ) 


where  the  Prandtl  number  is  Pr  ■  UCp/<  and  the  particle  Reynolds 
number  is  defined  as 


Re 
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The  thermal  conductivity  k  is  determined  from  Sutherland's  law, 
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where  S2  =  194  K,  T0  -  300  K,  and  <Q  =  2.6273  W/(m-K).  The  mean 
and  turbulent  heat  flux  vectors  may  be  written  as 


(A-34) 
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and 
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(A-35) 


where  the  turbulent  thermal  conductivity  is  defined  from 
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with  Preff  =  0.9.  The  mean  temperature  at  the  interface  between  the 
phases,  Tj , is  assumed  to  be  the  average  =  0.5  (T  +  TpS). 

The  Nobel-Abel  equation  of  state  for  the  gas  is 


pR„T 

pli-pv)  •  =  fizr  <*-”> 

m 

where  is  the  universal  gas  constant,  Wm  is  the  gas  molecular 
weight  and  n  is  the  covolume  factor  which  provides  a  correction  to  the 
perfect  gas  equation  of  state  for  large  density.  Other  required 
thermodynamic  properties  are 
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The  required  boundary  condition  for  Equation  (A-10)  at  the 
propellant  grain  surface  is  the  surface  heat  flux  relation, 


dJ 
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where  <q>  is  given  by  Equation  (A-29),  kp  is  the  particle  thermal 
conductivity,  and  $  is  the  heat  feedback  from  the  flame  identified  by 
Gough  (Reference  17)  as 


<d> 


^  Tps  V 


(A-42) 


where  dp  is  the  thermal  diffusivity  of  the  particles  and  Tpo  is 
taken  as  the  temperature  at  the  center  of  the  particle.  Following 
Reference  17,  the  heat  feedback  is  utilized  only  for  the  transient 
combustion  analysis  and  is  excluded  when  steady  state  combustion  is 
assumed. 
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